#include "Astronomy/SunMoonEphemeris.hpp"
#include "Astronomy/Celestial/Earth.hpp"
#include "Astronomy/Celestial/Moon.hpp"
#include "Astronomy/Celestial/Sun.hpp"
#include "Astronomy/Transformation.hpp"
#include "Astronomy/Constant.hpp"
#include <cmath>

double SunMoonEphemeris::ApparentEclipticLongitudeDiff(double T, int moonCount, int sunCount) noexcept
{
    double moonLongitude          = CelestialMoon::EclipticLongitude(T, moonCount);
    double moonVariationLongitude = CelestialMoon::VariationEclipticLongitude(T);

    double sunLongitude           = CelestialSun::GeometricEclipticLongitude(T, sunCount);
    double sunVariationLongitude  = CelestialSun::VariationEclipticLongitude(T);

    return moonLongitude + moonVariationLongitude - (sunLongitude + sunVariationLongitude + Constants::PI);
    //return moonLongitude + moonVariationLongitude - (sunLongitude + sunVariationLongitude);
}
double SunMoonEphemeris::GetEclipticLongitudeTime(double W) noexcept
{
    double v = 7771.37714500204;
    double t = (W + 1.08472) / v;
    t += (W - ApparentEclipticLongitudeDiff(t, 3, 3)) / v;
    v = CelestialMoon::Velocity(t) - CelestialEarth::Velocity(t); // v的精度0.5%，详见原文;
    t += (W - ApparentEclipticLongitudeDiff(t, 20, 10)) / v;
    t += (W - ApparentEclipticLongitudeDiff(t, -1, 60)) / v;
    return t;
}
double SunMoonEphemeris::GetEclipticLongitudeTimeLow(double W) noexcept
{
    double v  = 7771.37714500204;
    double t  = (W + 1.08472) / v;
    double t2 = t * t;
    t -= (-0.00003309 * t2 + 0.10976 * cos(0.784758 + 8328.6914246 * t + 0.000152292 * t2) + 0.02224 * cos(0.18740 + 7214.0628654 * t - 0.00021848 * t2) -
          0.03342 * cos(4.669257 + 628.307585 * t)) / v;

    double L = CelestialMoon::EclipticLongitude(t, 20) -
               (4.8950632 + 628.3319653318 * t + 0.000005297 * t * t + 0.0334166 * cos(4.669257 + 628.307585 * t) +
                0.0002061 * cos(2.67823 + 628.307585 * t) * t + 0.000349 * cos(4.6261 + 1256.61517 * t) - Transformation::ArcsecToRadians(20.5));
    v = 7771.38 - 914.0 * sin(0.7848 + 8328.691425 * t + 0.0001523 * t * t) - 179 * sin(2.543 + 15542.7543 * t) - 160 * sin(0.1874 + 7214.0629 * t);
    t += (W - L) / v;
    return t;
}